Hui-Walter Model

The multinomial distribution

Binomial (always with two possible outcomes):


Multinomial with two possible outcomes:


Multinomial with four possible outcomes:

Model Specification

model{
  Tally ~ dmulti(prob, N)
  
  # Test1- Test2-
    prob[1] <- (prev * ((1-se[1])*(1-se[2]))) + ((1-prev) * ((sp[1])*(sp[2])))

  # Test1+ Test2-
    prob[2] <- (prev * ((se[1])*(1-se[2]))) + ((1-prev) * ((1-sp[1])*(sp[2])))

  # Test1- Test2+
    prob[3] <- (prev * ((1-se[1])*(se[2]))) + ((1-prev) * ((sp[1])*(1-sp[2])))
  
  # Test1+ Test2+
    prob[4] <- (prev * ((se[1])*(se[2]))) + ((1-prev) * ((1-sp[1])*(1-sp[2])))

 

  prev ~ dbeta(1, 1)
  se[1] ~ dbeta(1, 1)
  sp[1] ~ dbeta(1, 1)
  se[2] ~ dbeta(1, 1)
  sp[2] ~ dbeta(1, 1)

  #data# Tally, N
  #monitor# prev, prob, se, sp, deviance
  #inits# prev, se, sp
}

twoXtwo <- matrix(c(48, 12, 4, 36), ncol=2, nrow=2)
twoXtwo
##      [,1] [,2]
## [1,]   48    4
## [2,]   12   36
library('runjags')

Tally <- as.numeric(twoXtwo)
N <- sum(Tally)

prev <- list(chain1=0.05, chain2=0.95)
se <- list(chain1=c(0.01,0.99), chain2=c(0.99,0.01))
sp <- list(chain1=c(0.01,0.99), chain2=c(0.99,0.01))

results <- run.jags('basic_hw.txt', n.chains=2)

[Remember to check convergence and effective sample size!]


results

Lower95 Median Upper95 SSeff psrf
prev 0.328 0.499 0.668 4250 2.288
prob[1] 0.367 0.462 0.558 13858 1.000
prob[2] 0.072 0.132 0.202 14200 1.000
prob[3] 0.018 0.055 0.104 9602 1.000
prob[4] 0.255 0.343 0.438 13555 1.000
se[1] 0.028 0.570 1.000 4564 15.155
se[2] 0.000 0.385 0.965 4585 13.619
sp[1] 0.000 0.461 0.970 4527 15.261
sp[2] 0.036 0.581 1.000 4593 13.641
deviance 12.304 15.134 21.447 8838 1.000






Label Switching

How to interpret a test with Se=0% and Sp=0%?

  • The test is perfect - we are just holding it upside down…

We can force se+sp >= 1:

  se[1] ~ dbeta(1, 1)
  sp[1] ~ dbeta(1, 1)T(1-se[1], )

Or:

  se[1] ~ dbeta(1, 1)T(1-sp[1], )
  sp[1] ~ dbeta(1, 1)

This allows the test to be useless, but not worse than useless.


Alternatively we can have the weakly informative priors:

  se[1] ~ dbeta(2, 1)
  sp[1] ~ dbeta(2, 1)

To give the model some information that we expect the test characteristics to be closer to 100% than 0%.

Or we can use stronger priors for one or both tests.

Priors

A quick way to see the distribution of a prior:

curve(dbeta(x, 1, 1), from=0, to=1)

qbeta(c(0.025,0.975), shape1=1, shape2=1)
## [1] 0.025 0.975

This was minimally informative, but how does that compare to a weakly informative prior for e.g. sensitivity?

curve(dbeta(x, 2, 1), from=0, to=1)

qbeta(c(0.025,0.975), shape1=2, shape2=1)
## [1] 0.1581139 0.9874209

qbeta(c(0.025,0.975), shape1=2, shape2=1)
## [1] 0.1581139 0.9874209

Or more accurately:

library("TeachingDemos")
hpd(qbeta, shape1=2, shape2=1)
## [1] 0.2236068 1.0000000

Credible vs confidence intervals:

  • For MCMC these are usually calculated using highest posterior density (HPD) intervals
  • Therefore there is a difference between:
    • qbeta(c(0.025,0.975), ...)
    • hpd(qbeta, ...)
  • Technically HPD intervals are credible intervals…

What about a more informative prior?

curve(dbeta(x, 20, 2), from=0, to=1)

qbeta(c(0.025,0.975), shape1=20, shape2=2)
## [1] 0.7618401 0.9882507
hpd(qbeta, shape1=20, shape2=2)
## [1] 0.7919691 0.9973994

Choosing a prior

What we want is e.g. Beta(20,1)

But typically we have median and 95% confidence intervals from a paper, e.g.:

“The median (95% CI) estimates of the sensitivity and specificity of the shiny new test were 94% (92-96%) and 99% (97-100%) respectively”

How can we generate a Beta( , ) prior from this?

The PriorGen package

“The median (95% CI) estimates of the sensitivity and specificity of the shiny new test were 94% (92-96%) and 99% (97-100%)”

library("PriorGen")
## Loading required package: rootSolve
findbeta(themedian = 0.94, percentile.value = 0.92)
## [1] "The desired Beta distribution that satisfies the specified conditions is: Beta( 429.95 27.76 )"
## [1] "Here is a plot of the specified distribution."
## [1] "Descriptive statistics for this distribution are:"
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.8906  0.9322  0.9401  0.9395  0.9474  0.9721 
## [1] "Verification: The percentile value 0.92 corresponds to the 0.05 th percentile"
hpd(qbeta, shape1=429.95, shape2=27.76)
## [1] 0.917172 0.960435

Note: themedian could also be themean


curve(dbeta(x, shape1=429.95, shape2=27.76))

Initial values

Part of the problem before was also that we were specifying extreme initial values:

se <- list(chain1=c(0.01,0.99), chain2=c(0.99,0.01))
sp <- list(chain1=c(0.01,0.99), chain2=c(0.99,0.01))

Let’s change these to:

se <- list(chain1=c(0.5,0.99), chain2=c(0.99,0.5))
sp <- list(chain1=c(0.5,0.99), chain2=c(0.99,0.5))

Analysing simulated data

This is useful to check that we can recover parameter values!

# Set a random seed so that the data are reproducible:
set.seed(2022-09-12)

sensitivity <- c(0.9, 0.6)
specificity <- c(0.95, 0.9)
N <- 1000
prevalence <- 0.5

data <- tibble(Status = rbinom(N, 1, prevalence)) %>%
  mutate(Test1 = rbinom(N, 1, sensitivity[1]*Status + (1-specificity[1])*(1-Status))) %>%
  mutate(Test2 = rbinom(N, 1, sensitivity[2]*Status + (1-specificity[2])*(1-Status)))

twoXtwo <- with(data, table(Test1, Test2))
Tally <- as.numeric(twoXtwo)

We know that e.g. the first test has Sensitivity of 90% and Specificity of 95% - so the model should be able to tell us that…

Practical Session 2

Points to consider

  1. What is the typical autocorrelation (and therefore effective sample size) of Hui-Walter models compared to the simpler models we were running earlier? Is there any practical consequence of this?

  2. How does changing the prior distributions for the se and sp of one test affect the inference for the other test parameters?

Summary

  • Hui-Walter models are seemingly magical, but:
    • They typically exhibit high autocorrelation
    • They may not converge, particularly with 1 population (see later!)
    • Need a larger sample for the same effective sample size
  • More informative priors for one test will
    • Improve identifiability of the model
    • Affect the posterior inference for the other test!